Baltimore releases detailed data on issues relevant to the city at https://data.baltimorecity.gov.

Detailed data on crime is released with a lag of just a week or so.

This allows us to get a good idea of the distribution of shootings in Baltimore.

This presentation uses R, a statistical programming language.

These slides are autogenerated during the data analytic process using ioslides, allowing for complete reproducibility.

All the code is embedded in this file but it won’t be displayed in every slide. For the full code you can go to https://github.com/peterphalen/ceasefire/mirecc-presentation, but check out https://peterphalen.com/ceasefire for a reader-friendly walkthrough.

The data on shootings comes from https://data.baltimorecity.gov/Public-Safety/BPD-Part-1-Victim-Based-Crime-Data/wsfq-mvij/data

url <- "https://raw.githubusercontent.com/peterphalen/ceasefire/
        master/BPD_Part_1_Victim_Based_Crime_Data.csv"
bpd <- readr::read_csv(url)
names(bpd)

##  [1] "CrimeDate"       "CrimeTime"       "CrimeCode"      
##  [4] "Location"        "Description"     "Inside/Outside" 
##  [7] "Weapon"          "Post"            "District"       
## [10] "Neighborhood"    "Longitude"       "Latitude"       
## [13] "Location 1"      "Premise"         "crimeCaseNumber"
## [16] "Total Incidents"

The dataset documents many types of crimes. We subset to shootings, i.e., non-lethal shootings and homicides with a firearm.

bpd <- subset(bpd, Description == "SHOOTING" | (Description == 
    "HOMICIDE" & Weapon == "FIREARM"))

There are usually many shootings per day, so we collapse the data down to daily counts.

bpd$CrimeDate <- as.Date(bpd$CrimeDate, format = "%m/%d/%Y")
daily <- bpd %>% group_by(CrimeDate) %>% summarise(shootings = n())
# fill missing dates, because many had no shootings
full.ts <- data.frame(CrimeDate = seq(daily$CrimeDate[1], 
    daily$CrimeDate[nrow(daily)], by = "day"))
daily <- full_join(full.ts, daily)
daily <- daily %>% group_by(CrimeDate) %>% mutate_all(funs(ifelse(is.na(.), 
    0, .)))

We need to represent weekly and yearly seasonality, in addition to the overall time trend.

library(lubridate)

daily$weekday <- weekdays(daily$CrimeDate)
daily$weekday <- factor(daily$weekday, levels = c("Monday", 
    "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday", 
    "Sunday"))

daily$day.of.year <- yday(daily$CrimeDate)

# the julian calendar is a simple system for
# numeric dates
daily$jul <- julian(daily$CrimeDate)

Poisson link function (outcome is count and mean is ~identical to sd)

Predictors: spline time trend, cyclical spline for yearly seasonality, and random effect for day of the week

library(rstanarm)
m1 <- stan_gamm4(shootings ~ 
            s(jul) +
            s(day.of.year, 
             bs="cc"), #cyclical constraint 
           random= ~ (1 | weekday),
           data=daily, cores=4,
           iter=500, family=poisson)

We draw the model predictions over the data in order to check the fit.

The time series is the sum of these three components:

Baltimore Ceasefire 365 is a city-wide call asking Baltimore residents to avoid having any murders through quarterly Ceasefires and Peace Challenges.

Ceasefires have been called four times per year since August 2017. These are “ceasefire weekends” but their impact often extends well beyond a few days (one lasted twelve).

This is a locally organized, non-governmental coalition. It is different from “Ceasefires” happening in most other cities.

Create a categorical variable for days during a Ceasefire.

# first day (Friday) of ceasefire weekends
ceasefire.initial <- as.Date(c("08/04/2017", "11/03/2017", 
    "02/02/2018", "05/11/2018", "08/03/2018", "11/02/2018", 
    "02/01/2019"), format = "%m/%d/%Y")

ceasefire.weekends <- lapply(ceasefire.initial, function(x) {
    seq(from = x, by = "day", length.out = 3)
})

ceasefire.weekends <- do.call("c", ceasefire.weekends)

# dummy variable
daily$ceasefire <- factor(ifelse(daily$CrimeDate %in% 
    ceasefire.weekends, 1, 0), labels = c("Regular Day", 
    "Ceasefire Weekend"))

m2 <- stan_gamm4(shootings ~ 
            ceasefire +
            s(jul) +
            s(day.of.year, 
             bs="cc"), #cyclical constraint 
           random= ~ (1 | weekday),
           data=daily, 
           cores=4,
           iter=500, 
           family=poisson)

plot(m2, regex_pars = "ceasefire") + 
scale_y_discrete(labels = "ceasefire effect") + 
xlim(c(-1.5, 0))

Model-estimated impact of ceasefire on Saturdays during the summer of 2018: